home *** CD-ROM | disk | FTP | other *** search
/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48_1 / leastsq < prev    next >
Internet Message Format  |  1995-03-31  |  13KB
  1. Path: seq!spell
  2. From: A.K. Donno <dnnant01@ucthpx.uct.ac.za>
  3. Subject:  v02i003:  leastsq - Least squares curve fitting, Part01/01
  4. Newsgroups: comp.sources.hp48
  5. Keywords: Least squares
  6. Organization: University of Cape Town
  7. Followup-To: comp.sys.hp48
  8. Summary: A program that will fit a polynomial to a set a data points
  9. Approved: spell@seq.uncwil.edu
  10.  
  11. Checksum: 2430449214 (verify with brik -cv)
  12. Submitted-by: A.K. Donno <dnnant01@ucthpx.uct.ac.za>
  13. Posting-number: Volume 2, Issue 3
  14. Archive-name: leastsq/part01
  15.  
  16. BEGIN_DOC leastsq.doc
  17. The feature that I use most in the HP48 stats menu is that of curve fitting,
  18. however if you have data that does not fit the standard curve types offered
  19. by the calculator, a good fit cannot be obtained.
  20.  
  21. The following program fits a polynomial equation to a set of data points in a
  22. similar fashion to the stats menu. It takes a matrix off the stack that has at
  23. least two columns in it (this is not checked), allows you to specify which
  24. column is x data and which is y data, and then fits a curve to the data using
  25. a polynomial of an order that you specify. Please note that more than 2
  26. columns are allowed, but the program only works with the two that are
  27. selected.
  28.  
  29. The program uses a standard least squares fit, but operates fairly quickly
  30. (4.5 mins to fit a 15th order to 21 data pairs). Note that an option allows
  31. you to specify whether the program must plot the fitted curve, or just return
  32. the equation. The program does not recalculate the polynomial, unless you
  33. change the x or y columns or the order of the polynomial required.
  34.  
  35. The following menu is used:
  36.  
  37. XCOL   - takes an integer off the stack and makes that column the x data
  38.          column
  39. YCOL   - takes an integer off the stack and makes that column the y data
  40.          column
  41. ORDER  - takes an integer off the stack and makes it the order of the poly-
  42.          nomial that is required
  43. PLOT   - calculates (if necessary) and plots the fitted polynomial
  44. FIT    - calculates (if necessary) and returns the polynomial that fits the
  45.          data
  46. EXIT   - returns back to the LEAST directory CST menu (you can put whatever
  47.          you want in this menu although a suggested menu is given)
  48.  
  49. The following is a list of the programs that make up LEAST:
  50.  
  51. LEAST  - main program (menus ...)
  52. PLOT   - plotting function
  53. FIT    - curve fitter (brains)
  54. POLY   - converts a list of numbers on the stack into a algebraic polynomial
  55.          of the order given in level one (the numbers are the coefficients)
  56. SETRNG - sets the range for the plot
  57. XYR    - redraws the top line of the LEAST screen (XCOL: 1 YCOL: 2 ORDER: 3)
  58. DELEQ  - deletes the current equation EQ
  59.  
  60. The following variables are produced by LEAST:
  61.  
  62. ORDR  - order of the polynomial
  63. YC    - y data column
  64. XC    - x data column
  65. dDAT  - summation data (same as the stats menu)
  66. EQ    - current equation
  67.  
  68. Checksum : # 63743d
  69. Bytes    : 1950.5
  70.  
  71. I hope that this program is of use to others, and would appreciatiate it if
  72. you would mail any comments to me.
  73. END_DOC
  74.  
  75. Here is the ASC file for the program followed by the source:
  76. (decode with ASC->)
  77.                                     
  78. BEGIN_ASC leastsq.asc
  79. %%HP: T(3)A(D)F(.);
  80. "69A20FF7F390000000303435453047A2084E2050C45414354584E204005C4F44
  81. 584E203064944547A20B213047A20C2A20F00000555257454D9D20E163247A20
  82. 84E2020541584E20400505142584E2040F425442584E2020953484E202085348
  83. 4E204058441445B2130EFE0293632B2130B213047A20C2A20D000054859445D9
  84. D20E1632041A193632B2130B2130B2130F0100504454C4541550D9D20E16323C
  85. E224563284E202054159763282EC1683A2D9AE1AFE22D9D204563284E2020541
  86. 597632EFE02B21305DF2293632B2130970003085952530D9D20E1632916C11C4
  87. 32D6E205066C6167637E16323392010000000000005495D2C133920100000000
  88. 00006495D2C13392010000000000007495D2C13392010000000000008495D2C1
  89. 3392010000000000009495D2C13392010000000000000595D2C1C2A20F000085
  90. 36F6C6A384E20208534B0BC176BA1C2A2011000029536F6C6A376BA184E20209
  91. 534B0BC176BA1C2A203100002F427465627A376BA184E2040F4254425B0BC176
  92. BA1C2A2033000A0F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F576BA1
  93. 9C2A2485A19C2A24A5A1D6E205066C6167637F76C1EF53293632B21305D10060
  94. 35544525E47460D9D20E1632E7EF184E202085346C7D1E7EF184E202095346C7
  95. D11C432D6E204087D61687D6E204097D61687E16324BEF184E202085346C7D14
  96. BEF184E202095346C7D11C432D6E204087D696E6D6E204097D696E6E1632D6E2
  97. 04087D61687D6E204087D696E6126E1D6E204097D61687D6E204097D696E6146
  98. E1EF532EF53293632B2130811004005F4C49540D9D20E16321C432D6E201046E
  99. 16320B3A25D2C14563284E20108597632D6E201046D20B1EEDA1D6E2010469C2
  100. A290DA14B2A20A132D6E201046DBBF14563284E20108597632D6E201046D20B1
  101. EEDA176BA1683A208332EF53293632B21308C0003064944530D9D20E163284E2
  102. 040584414458B9C1B7FC18DBF18DBF11C432D6E201046E16323CE22D6E201046
  103. 84E2040F4254425D5CE1AFE22D9D20D6E20104684E2040F42544259C2A276BA1
  104. ED2A2387C19C2A2681D1ED2A284E2040F42544259C2A276BA10A132D6E201096
  105. 9C2A2D6E2010460A132D6E2010A6D6E2010A6D6E201096ED2A2387C184E20405
  106. 8441445D6E2010A684E20208534ED2A2387C16C7D1D6E2010969C2A290DA1D20
  107. B1704D1C4232C42329C2A2D6E2010460A132D6E20109684E204058441445D6E2
  108. 0109684E20209534ED2A2387C16C7D1C4232D6E2010469C2A2ED2A2387C1900D
  109. 1DBBF178BF178BF1293D1DBBF1EEDA1872B1DBBF1293D1EEDA1DBBF1EEDA1B7F
  110. C18DBF184E2040F425442584E204005F4C495B21305BF22C2A20D5000542525F
  111. 425A302F42746562702E30205F696E64737A0F5F5F5F5F5F5F5F5F5F5F5F5F5F
  112. 5F5F5F5F5F5F5F5F55DF224563284E2020541597632DCC02EF53293632B21307
  113. D2004005C4F44540D9D20E1632F52E184E2040584414458B9C1B7FC18DBF18DB
  114. F11C432D6E201046E16329C2A2D6E2010460A132D6E20109684E204058441445
  115. D6E20109684E20208534ED2A2387C16C7D184E204058441445D6E20109684E20
  116. 209534ED2A2387C16C7D1E97C1A92E1A13E1C42326C1E1091E1EF53293632B21
  117. 30EF00050C45414354550D9D20E163284E20504454C45415B2DF116DF1858A19
  118. C2A24563284E2020853497632DCC02ED2A24563284E2020953497632DCC023F2
  119. A24563284E2040F425442597632DCC0284E203085952547A2047A20C2A20D000
  120. 08534F4C447A20D9D20E16324563284E2020853497632DCC0284E20504454C45
  121. 41584E203085952593632B2130D9D20E16324563284E2020853497632DCC0284
  122. E20504454C4541584E203085952593632B2130D9D20E163284E2020853484E20
  123. 3085952593632B2130B2130B213047A20C2A20D00009534F4C447A20D9D20E16
  124. 324563284E2020953497632DCC0284E20504454C4541584E203085952593632B
  125. 2130D9D20E16324563284E2020953497632DCC0284E20504454C4541584E2030
  126. 85952593632B2130D9D20E163284E2020953484E203085952593632B2130B213
  127. 0B213047A20C2A20F0000F42544542547A20D9D20E16324563284E2040F42544
  128. 2597632DCC0284E20504454C4541584E203085952593632B2130D9D20E163245
  129. 63284E2040F425442597632DCC0284E20504454C4541584E203085952593632B
  130. 2130D9D20E163284E2040F425442584E203085952593632B2130B2130B213047
  131. A20C2A20D000005C4F445D9D20E16323CE224563284E202054159763282EC168
  132. 3A2279E1AFE22D9D2084E206035544525E47484E2030649445B21305DF223CE2
  133. 24563284E202054159763282EC1ED2A2279E1AFE22D9D2084E202054159C2A24
  134. 85A1ED2A2F17A1B21305BF22D9D20166E184E204005C4F445AB2E1B21305DF22
  135. 84E203085952593632B2130B213047A20C2A20B0000649445D9D20E16323CE22
  136. 4563284E202054159763282EC1683A2279E1AFE2284E20306494455DF223CE22
  137. 4563284E202054159763282EC1ED2A2279E1AFE22D9D2084E202054159C2A248
  138. 5A1ED2A2F17A1B21305BF22D9D204563284E202054159763204B02B21305DF22
  139. 84E203085952593632B2130B213047A20C2A20D000054859445D9D20E16324B2
  140. A26911293632B2130B2130B2130D511293632B2130FF8F"
  141. END_ASC
  142.  
  143.  
  144. BYTES: #F8FFh 1952.5
  145.  
  146. BEGIN_UU leastsq.uue
  147. begin 644 leastsq.bin
  148. M2%!(4#0X+466*O!_/PD````#0U-4`W0J@.0"!4Q%05-42"Y``,7T1(7D`@-&\
  149. M251T*K`2`W0JP*("#P``525U5-39`AXV0J<"2"X@4!2%Y`($4%!!4D@N0/`D"
  150. M122%Y`("64-(+B"`-83D`@2%1$%4*S'@[R`Y-K(2`RLQ0*<"+"K0``!%6$E4P
  151. MG2W@82-`H9%C(RLQL!(#*S'P$``%1$5,15$%G2W@82/#+D)E(T@N(%`4E6<CO
  152. M*,YA."J=ZJ'O(ITM0&4C2"X@4!259R/^#K(2`]4ODF,C*S&0!P`#6%E2`YTMT
  153. MX&$C&<813"-M+E!@QA9V-N=A(S,I$````````$59+1PS*1````````!&62T<I
  154. M,RD0````````1UDM'#,I$````````$A9+1PS*1````````!)62T<,RD0````8
  155. M````4%DM'"PJ\```6&-O;#I(+B"`-;2P'&>KP:("$0``DC7VQJ9SMAI(+B"0`
  156. M-;2P'&>KP:("$P``\B1'5B:G<[8:2"Y`\"1%)+6P'&>KP:(",P"@\/7U]?7U!
  157. M]?7U]?7U]?7U]?7U]?7U]76V&LFB0E@:R:)"6AIM+E!@QA9V-O=G'/XUDF,CC
  158. M*S%0'0`&4T544DY'!ITMX&$C?OZ!Y`("6$/&U^'G'T@N()`U9'P=P332Y@($M
  159. M>&UA>&TN0)#7%H;G82.T_H'D`@)80\;70>L?2"X@D#5D?!W!--+F`@1X;6ENT
  160. M;2Y`D->6YN9A(VTN0(#7%H;7Y@($>&UI;B'FT>8"!'EM87AM+D"0UY;F%F0>Z
  161. M_C7B7R,Y-K(2`Q@!0`#UQ)1%T-D"'C823"-M+A!`YF$CL*-2+1Q4-H+D`@%8Y
  162. M>3;2Y@(!9"VPX=X:;2X00)8L*@FM02LJH#'2Y@(!9+W[064C2"X0@)5G(VTNI
  163. M$$#6`AONK7&V&H:C`C@C_C628R,K,8`,``-&250#G2W@82-(+D!02!1$A9L<O
  164. M>\^!O1_8^Q%,(VTN$$#F82/#+M+F`@%D2"Y`\"1%)-7%'OHNTMD";2X00(;DA
  165. M`@1/4D12R:)RMAK>HC)X',FB8A@=WJ*"Y`($3U)$4LFB<K8:H#'2Y@(!:<FB"
  166. MTN8"`62@,=+F`@%J;2X0H-;F`@%IWJ(R>!Q(+D!02!1$U>8"`6I(+B"`->0M,
  167. M*H/'87P=;2X0D)8L*@FMT0(;!]3!)"-,,I(L*FTN$$`&&B-M+A"0AN0"!(5$P
  168. M051M+A"0AN0"`EE#WJ(R>!S&U\$D(VTN$$"6+"K>HC)X'`G0T;L?A_MQN!^22
  169. MT]&['^ZM@2<;O?LA.1WNK=&['^ZML?<<V/N!Y`($3U)$4D@N0`#UQ)2U$@.UB
  170. M+\*B`ET`4"0E]22E`_(D1U8F!^(#`O66YD8WI_#U]?7U]?7U]?7U]?7U]?7U%
  171. M]?7U]?55_2)4-H+D`@)%47DVTLP@_C628R,K,7`M``103$]4!)TMX&$C7^*!W
  172. MY`($A41!5+C)L?<<V/N!O1_!--+F`@%D'C:2+"IM+A!`!AHC;2X0D(;D`@2%-
  173. M1$%4;2X0D(;D`@)80]ZB,G@<QM>!Y`($A41!5&TN$)"&Y`("64/>HC)X',;7M
  174. MX7D<FN*A,1Y,,F(<'I#AX5\C.3:R$@/^`%#`5!0T1570V0(>-H+D`@5$14Q%W
  175. M42O]$=8?6*B1+"I4-H+D`@)80WDVTLP@WJ)"92-(+B"0-91G(\T,,B\J5#:"2
  176. MY`($3U)$4GDVTLP@2"XP@)4E1:<"="K`H@(-`(`U],1$IP*=+>!A(U0V@N0"V
  177. M`EA#>3;2S"!(+E!`5,14%(7D`@-865(Y-K(2`YTMX&$C5#:"Y`("6$-Y-M+,.
  178. M($@N4$!4Q%04A>0"`UA94CDVLA(#G2W@82-(+B"`-83D`@-865(Y-K(2`RLQV
  179. ML!(#="K`H@(-`)`U],1$IP*=+>!A(U0V@N0"`EE#>3;2S"!(+E!`5,14%(7D1
  180. M`@-865(Y-K(2`YTMX&$C5#:"Y`("64-Y-M+,($@N4$!4Q%04A>0"`UA94CDVL
  181. MLA(#G2W@82-(+B"0-83D`@-865(Y-K(2`RLQL!(#="K`H@(/`/`D150D1:<"+
  182. MG2W@82-4-H+D`@1/4D12>3;2S"!(+E!`5,14%(7D`@-865(Y-K(2`YTMX&$CW
  183. M5#:"Y`($3U)$4GDVTLP@2"Y00%3$5!2%Y`(#6%E2.3:R$@.=+>!A(T@N0/`DS
  184. M122%Y`(#6%E2.3:R$@,K,;`2`W0JP*("#0``Q?1$U=D"'C8R[")4-H+D`@)%6
  185. M47DV@N(<AJ,BEQ[Z+M+9`D@N8#!51"7E=(3D`@-&250K,5#](L,N0F4C2"X@B
  186. M4!259R,HSN$M*G+IH>\BG2V`Y`("15')HD)8&MZB\G$:*S%0^R*=+1!F'D@N@
  187. M0`#%]$2E*QXK,5#](D@N,("5)95C(RLQL!(#="K`H@(+`&"41-79`AXV,NPB'
  188. M5#:"Y`("15%Y-H+B'(:C(I<>^BZ"Y`(#1DE4U2\R[")4-H+D`@)%47DV@N(<:
  189. MWJ(BEQ[Z+M+9`D@N(%`4E2PJA*7A+2H?I[$2`[4OTMD"5#:"Y`("15%Y-@*TB
  190. M("LQ4/TB2"XP@)4EE6,C*S&P$@-T*L"B`@T`4(251-79`AXV0BLJEA&28R,KU
  191. .,;`2`RLQT!4A.3:R$@.R#
  192. ``
  193. end
  194. END_UU
  195.  
  196. BEGIN_RPL leastsq.rpl
  197. %%HP: T(3)A(D)F(.);
  198. DIR
  199.   LEAST
  200.     \<< DELEQ CL\GS \GS+
  201. CLLCD 1 'XC' STO 2
  202. 'YC' STO 3 'ORDR'
  203. STO XYR { { "XCOL"
  204. {
  205.       \<< 'XC' STO
  206. DELEQ XYR
  207.       \>>
  208.       \<< 'XC' STO
  209. DELEQ XYR
  210.       \>>
  211.       \<< XC XYR
  212.       \>> } } {
  213. "YCOL" {
  214.       \<< 'YC' STO
  215. DELEQ XYR
  216.       \>>
  217.       \<< 'YC' STO
  218. DELEQ XYR
  219.       \>>
  220.       \<< YC XYR
  221.       \>> } } {
  222. "ORDER" {
  223.       \<< 'ORDR' STO
  224. DELEQ XYR
  225.       \>>
  226.       \<< 'ORDR' STO
  227. DELEQ XYR
  228.       \>>
  229.       \<< ORDR XYR
  230.       \>> } } {
  231. "PLOT"
  232.       \<<
  233.         IF 'EQ'
  234. VTYPE -1 ==
  235.         THEN SETRNG
  236. FIT
  237.         END
  238.         IF 'EQ'
  239. VTYPE 2 ==
  240.         THEN EQ 1
  241. DISP 2 WAIT
  242.         ELSE
  243. FUNCTION PLOT GRAPH
  244.         END XYR
  245.       \>> } { "FIT"
  246.       \<<
  247.         IF 'EQ'
  248. VTYPE -1 ==
  249.         THEN FIT
  250.         END
  251.         IF 'EQ'
  252. VTYPE 2 ==
  253.         THEN EQ 1
  254. DISP 2 WAIT
  255.         ELSE 'EQ'
  256. RCL
  257.         END XYR
  258.       \>> } { "EXIT"
  259.       \<< 0 MENU
  260.       \>> } } TMENU
  261.     \>>
  262.   PLOT
  263.     \<< ERASE \GSDAT
  264. SIZE OBJ\-> DROP DROP
  265. \-> d
  266.       \<< 1 d
  267.         FOR i \GSDAT
  268. i XC 2 \->LIST GET
  269. \GSDAT i YC 2 \->LIST
  270. GET R\->C C\->PX PIXON
  271.         NEXT DRAX
  272. DRAW
  273.       \>>
  274.     \>>
  275.   FIT
  276.     \<< \GSDAT SIZE
  277. OBJ\-> DROP DROP \-> d
  278.       \<<
  279.         IF d ORDR >
  280.         THEN d ORDR
  281. 1 + 2 \->LIST 1 CON 2
  282. ORDR 1 +
  283.           FOR i 1 d
  284.             FOR j j
  285. i 2 \->LIST \GSDAT j XC
  286. 2 \->LIST GET i 1 - ^
  287. PUT
  288.             NEXT
  289.           NEXT 1 d
  290.           FOR i
  291. \GSDAT i YC 2 \->LIST
  292. GET
  293.           NEXT d 1
  294. 2 \->LIST \->ARRY SWAP
  295. DUP DUP TRN SWAP *
  296. INV SWAP TRN * SWAP
  297. * OBJ\-> DROP ORDR
  298. POLY
  299.         ELSE
  300. "ERROR: Order > Points
  301. ______________________"
  302.         END 'EQ'
  303. STO
  304.       \>>
  305.     \>>
  306.   POLY
  307.     \<< \-> d
  308.       \<< -3 CF 'X' d
  309. ^ * d 1 - 0
  310.         FOR d SWAP
  311. 'X' d ^ * + -1
  312.         STEP
  313.       \>>
  314.     \>>
  315.   SETRNG
  316.     \<< MAX\GS XC GET
  317. MAX\GS YC GET \-> xmax
  318. ymax
  319.       \<< MIN\GS XC GET
  320. MIN\GS YC GET \-> xmin
  321. ymin
  322.         \<< xmax xmin
  323. XRNG ymax ymin YRNG
  324.         \>>
  325.       \>>
  326.     \>>
  327.   XYR
  328.     \<< RCLF \-> flags
  329.       \<< -45 CF -46
  330. CF -47 CF -48 CF
  331. -49 CF -50 CF
  332. "Xcol:" XC \->STR +
  333. " Ycol:" + YC \->STR
  334. + " Order:" + ORDR
  335. \->STR +
  336. "
  337. ______________________"
  338. + 1 DISP 1 FREEZE
  339. flags STOF
  340.       \>>
  341.     \>>
  342.   DELEQ
  343.     \<<
  344.       IF 'EQ' VTYPE
  345. -1 \=/
  346.       THEN 'EQ'
  347. PURGE
  348.       END
  349.     \>>
  350.   CST { LEAST PLOT
  351. FIT { } { "PURGE"
  352.     \<< { EQ PPAR
  353. ORDR YC XC \GSDAT }
  354. PURGE
  355.     \>> } { "EXIT"
  356.     \<< HOME
  357.     \>> } }
  358. END
  359. END_RPL
  360. ==============================================================================
  361. Anthony Donno                          |
  362. Department of Electrical Engineering   |   Caesar entered on his head
  363. University of Cape Town                |   he had a helmet on each foot
  364. New South Africa                       |   he had a sandal in his hand
  365.                                        |   he had his trusty sword to boot
  366. e-mail : dnnant01@ucthpx.uct.ac.za     |
  367. phone  : +27 +21 6892420               |
  368. ==============================================================================
  369.  
  370.  
  371.